En vez de evaluarla analíticamente, podemos simular de esta distribución utilizando las muestras de \(\theta\) que obtuvimos con MCMC.
Para cada observación \(y_i\) (con \(i \in \{1, ..., N\}\)), y con cada muestra de la posterior \(\theta_s\) (con \(s \in \{1, ..., S\}\)), se simula al menos una observación de la variable respuesta (\(\tilde{y}_i^s\)).
Si queremos utilizar esta distribución para verificar el modelo, tenemos que simular nuevos datos utilizando los mismos valores de las predictoras (\(x\)) asociados a cada observación (\(x_i\)).
Modelo de incendios
Ejemplo con el modelo de la cantidad de incendios por año
Code
# Cargamos datosdatos <-read.csv(here::here("R_curso_modelos", "datos","barbera_data_fire_total_climate.csv"))# Compilamos el modelomodel <-cmdstan_model(here::here("R_curso_modelos", "modelos", "nfuegos.stan"))# Creamos una lista nombrada para pasarle los datos. Los nombres de cada# elemento de la lista tienen que ser exactamente los que definimos en # la sección data {}stan_data <-list(N =nrow(datos),y = datos$fires,x = datos$fwi)# Y muestreamos la posteriorfit_mcmc <- model$sample(data = stan_data,chains =4, parallel_chains =4,iter_warmup =1000,iter_sampling =1000)
Processing csv files: /tmp/RtmpVB3dXY/nfuegos-202505052112-1-8d555f.csv, /tmp/RtmpVB3dXY/nfuegos-202505052112-2-8d555f.csv, /tmp/RtmpVB3dXY/nfuegos-202505052112-3-8d555f.csv, /tmp/RtmpVB3dXY/nfuegos-202505052112-4-8d555f.csv
Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.
Checking sampler transitions for divergences.
No divergent transitions found.
Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.
Effective sample size satisfactory.
Split R-hat values satisfactory all parameters.
Processing complete, no problems detected.
Code
# Extraemos las muestras de alpha y beta, en formato data.framed <- fit_mcmc$draws(variables =c("alpha", "beta"), format ="draws_df")colnames(d) <-c("a", "b")# Ahora podemos simular desde la distribución predictiva posterior.# Para cada observación simularemos S = 4000 muestras de la predictiva posterior. # Guardaremos las simulaciones en una matriz de N * S:S <-nrow(d) # número de muestras de la posteriorN <-nrow(datos) # número de observacionesysim <-matrix(NA, N, S)# Simulamosfor (s in1:S) {# calculamos lambda para la muestra "s" lambda <-exp(d$a[s] + d$b[s] * datos$fwi)# es un vector de largo N, porque datos$fwi también lo es.# llenamos una columna entera en ysim, que equivale a un set de datos simulado ysim[, s] <-rpois(N, lambda)}# Ahora cada fila tiene muestras de la distribución predictiva posterior de y # para el valor correspondiente de fwi
Code
fila <-sample(1:N, 1)rr <-c(datos$fires[fila], ysim[fila, ]) |>range()hist(ysim[fila, ], xlim = rr, main =paste("Fila", fila),xlab ="Número de incendios", ylab ="Frecuencia")abline(v = datos$fires[fila], lty =2, col ="red", lwd =1.5)
Una vez obtenida la (muestra de la) distribución predictiva posterior, necesitamos evaluar si nuestros datos parecen salidos de ahí.
Para ello utilizamos la Transformación Integral de Probabilidad (PIT, de Probability Integral Transform):
sea \(Y\) una variable aleatoria continua y \(F_Y\) su función de probabilidad acumulada, la variable aleatoria \(Z\), definida como \(Z = F_Y(Y)\), sigue una distribución uniforme entre 0 y 1.
Esto implica que si nuestras observaciones se comportan como muestras de la distribución predictiva posterior, sus valores de probabilidad acumulada en conjunto siguen una distribución uniforme entre 0 y 1.
A diferencia de los distintos tipos de residuos que suelen calcularse, esto aplica sin importar qué características tenga la distribución predictiva posterior.
Sin embargo, esto no se cumple si \(Y\) es discreta, y en ese caso hay que aleatorizar los valores de probabilidad acumulada (\(Z\)) para que se distribuyan de forma uniforme.
Estos valores de probabilidad acumulada (aleatorizados en caso de variables discretas) son calculados por el paquete DHARMa (DHARMa residuals), pero también aparecen en INLA, loo y bayesplot bajo el nombre de PIT values.
Ejemplo en R
Code
N <-500y <-rnorm(N)p <-pnorm(y)hist(p, breaks =30)
Code
# N estadísticos de orden de una uniforme:q_unif <-ppoints(N)plot(sort(p) ~ q_unif, pch =19, col =rgb(0, 0, 0, 0.1),ylab ="Cuantiles observados", xlab ="Cuantiles esperados")abline(0, 1, col ="red")
Code
# ¿Y si las muestras son una mezcla proveniente de distintas distribuciones?y1 <-rnorm(N)y2 <-rgamma(N, 1)par(mfrow =c(1, 2))plot(density(y1))plot(density(y2, from =0))
q_unif <-ppoints(N)plot(sort(p) ~ q_unif, pch =19, col =rgb(0, 0, 0, 0.1),ylab ="Cuantiles observados", xlab ="Cuantiles esperados")abline(0, 1, col ="red")
Para transformar estos valores de probabilidad acumulada en una variable con distribución uniforme, construimos nuestro valor \(Z\) (PIT) de la siguiente manera: \[Z_i \sim \text{uniforme}(F_Y(y_i- 1), F_Y(y_i)).\]
Es decir, en vez de quedarnos con la probabilidad acumulada (\(F_Y(y_i)\)), tomamos una muestra de una distribución uniforme limitada entre ese valor (por arriba) y la probabilidad acumulada de \(y_i - 1\) (\(F_Y(y_i- 1)\)).
PIT para variable discreta
Code
y <-rpois(N, 10)freqs <-table(y) / Nvals <-as.numeric(names(freqs))barplot(freqs ~ vals)
Code
p <-ppois(y, lambda =10) # Probabilidad acumulada en ypl1 <-ppois(y-1, lambda =10) # Probabilidad acumulada en y-1pit <-runif(N, pl1, p)par(mfrow =c(1, 2))hist(p, breaks =20, main =NA, xlab ="Probabilidad acumulada")hist(pit, breaks =20, main =NA, xlab ="Probabilidad acumulada\naleatorizada")
Cuando evaluamos ajustes bayesianos, la distribución estimada por el modelo generalmente varía entre observaciones (sigue la misma estructura, e.g., poisson, pero tiene otros parámetros). En general aproximamos la acumulada mediante la función de probabilidad acumulada empírica (en R, la función ecdf).
Code
N <-nrow(datos)fila <-sample(1:N, 1)# Obtenemos la función de distribución acumulada empírica de la distribución # predictiva posterior para la fila elegidaecdf_post <-ecdf(ysim[fila, ])# al evaluarla, en el valor observado nos devuelve la probabilidad acumulada(prob <-ecdf_post(datos$fires[fila]))
[1] 0.14425
Code
# Calculamos el PIT:(pit <-runif(1, min =ecdf_post(datos$fires[fila] -1),max =ecdf_post(datos$fires[fila])))
Pero DHARMa calcula los PIT por nosotros; sólo necesitamos darle una matriz de datos simulados.
Code
res <-createDHARMa(simulatedResponse = ysim,observedResponse = datos$fires,integerResponse =TRUE# se da cuenta solo, pero por las dudas)plot(res)
Code
hist(res$scaledResiduals)
Code
# Podemos ver si hay alguna tendencia en función de una predictoraplotResiduals(res, form = datos$fwi)
El modelo admite menor variabilidad que la encontrada en los datos.
Este modelo debe ser mejorado, pero podemos usar los residuos para ver si nos falta agregar alguna variable importante:
Code
# Agregamos los PIT calculados por DHARMa a la tabla de datosdatos$pit <- res$scaledResidualsggplot(datos, aes(pp, pit)) +geom_point(size =1.5) +labs(y ="Prob. acumulada", # o "residuos dharma"x ="Precipitación (mm)") +geom_smooth(method ="gam")
Code
# No parece haber una tendencia, lo cual sugiere que no ganaríamos mucho # incluyendo la precipitación.
Esa es la forma más robusta de comparar la predictiva posterior con los datos, pero también podemos usar al paquete bayesplot, de nuestros amigos de Stan, para visualizar el ajuste de otras formas.
Code
# Trasponemos la matriz de datos simulados porque a bayesplot le gustan las # iteraciones de la posterior como filas (las teníamos en las columnas)ysim_t <-t(ysim)# posterior predictive intervals en función de algoppc_intervals(y = datos$fires, # observacionesyrep = ysim_t, # simulaciones de la predictiva posteriorx = datos$fwi, # alguna predictora de interésprob =0.9# probabilidad de los intervalos) +labs(x ="FWI", y ="Número de incendios")
Code
# Y la densidad marginal de la respuesta, según lo observado (y) y los sets de # datos simulados (yrep). Como graficará una curva por set de datos simulado, # será mejor elegir al azar unos pocos (menos que 4000)use <-sample(1:S, 100) # elige 100 al azarppc_dens_overlay(y = datos$fires, yrep = ysim_t[use, ])
Con estas herramientas podemos evaluar dónde falla el modelo y así intentar mejorarlo.
Cuánto hay que mejorarlo, depende del contexto.
Interpretación de modelos
Si el modelo se adecúa razonablemente a los datos, podemos comenzar a evaluar qué implica la posterior estimada. (Aunque a veces, visualizar el ajuste también puede ayudar a identificar potenciales problemas.)
En algunos casos los parámetros en sí resultan de interés, pero en otros casos, sólo nos interesan cantidades derivadas de ellos.
Ejemplo: nos interesa evaluar cómo varía \(\lambda\) en función del FWI.
Cada muestra \(s\) de la posterior (\(\{\alpha_s, \beta_s\}\)) define una función \[\lambda_s = f_s(\text{FWI}) = \exp(\alpha_s + \beta_s \ \text{FWI}).\]
Si nos interesa describir esa función considerando la incertidumbre en la estimación, podemos evaluarla sobre una secuencia de \(\text{FWI}\) utilizando todas o algunas muestras de la posterior.
Importante: si guardamos las muestras en una matriz, donde cada fila es una iteración, no podemos mezclar valores de distintas filas. Esto es porque estaríamos rompiendo la estructura de correlaciones de la posterior conjunta.
Code
nrep <-200# largo de la secuencia de FWIfwi_seq <-seq(min(datos$fwi), max(datos$fwi), length.out = nrep)# Guardaremos la función evaluada sobre la secuencia de fwi en una matriz, # donde cada columna corresponderá a una muestra de la posteriorlambda_mat <-matrix(NA, nrep, S)for (s in1:S) { lambda_mat[, s] <-exp(d$a[s] + d$b[s] * fwi_seq)}# Podemos graficar algunas curvas ("spaghetti plot").# Ordenaremos la matriz para graficarla con ggplotncurves <-50use <-sample(1:S, ncurves) # elegimos 150 curvas al azarlambda_mat_sub <-as.data.frame(lambda_mat[, use]) lambda_mat_sub$fwi <- fwi_seq# Elongamos la matrizllong <-pivot_longer( lambda_mat_sub, cols =all_of(1:ncurves), names_to ="sim", values_to ="lambda")ggplot(llong, aes(fwi, lambda, group = sim)) +geom_line(alpha =0.1) +labs(x ="FWI", y =expression(lambda))
Code
# Lambda es una cantidad derivada de alpha, beta y el fwi. Todo lo que se calcula# a partir de los parámetros tiene su propia distribución posterior. Y acá # tenemos, además, una distribución posterior de lambda para cada valor de # fwi que evaluemos.# Otra forma de graficar la función es resumir la distribución de cada lambda,# es decir, la distribución de lambda marginal a alpha y beta para cada x.pred <-data.frame(fwi = fwi_seq,lambda_mean =rowMeans(lambda_mat),lambda_lwr =apply(lambda_mat, 1, quantile, probs =0.025),lambda_upr =apply(lambda_mat, 1, quantile, probs =0.975))# apply() es una forma concisa de loopear. equivale a hacerfor (i in1:N) { pred$lambda_lwr[i] <-quantile(lambda_mat[i, ], probs =0.025)}# Es decir, aplica la función quantile, con argumento probs = 0.025, sobre el # array llamado "lambda_mat", iterando sobre la dimensión 1 (filas)# Ahora graficamos la predicción media y su incertidumbreggplot(pred, aes(fwi, lambda_mean, ymin = lambda_lwr, ymax = lambda_upr)) +geom_ribbon(alpha =0.4, color =NA, fill ="orange") +geom_line() +labs(x ="FWI", y =expression(lambda)) +# Agregamos los datosgeom_point(aes(fwi, fires), data = datos, inherit.aes = F,alpha =0.7, size =2)
Code
# La cinta muestra los percentiles 2.5 y 97.5 % de la distribución posterior de# lambda. Pero también podemos explorar los percentiles extremos de la# distribución predictiva posterior de y.# Para eso, simulamos nuevas observaciones sobre la secuencia de FWI.# (Antes, para el chequeo posterior, habíamos simulado nuevos valores de y# usando los valores observados de FWI, no sobre una secuencia.)ysim_seq <-matrix(NA, nrep, S)for (s in1:S) { ysim_seq[, s] <-rpois(nrep, lambda = lambda_mat[, s])}# Obtenemos los cuantiles extremos:pred$y_lwr <-apply(ysim_seq, 1, quantile, probs =0.025)pred$y_upr <-apply(ysim_seq, 1, quantile, probs =0.975)# Ahora graficamos la predicción media y su incertidumbreggplot(pred, aes(fwi, lambda_mean, ymin = lambda_lwr, ymax = lambda_upr)) +# Intervalo para ygeom_ribbon(aes(fwi, ymin = y_lwr, ymax = y_upr), inherit.aes = F, alpha =0.3, color =NA, fill ="blue") +# Intervalo para lambdageom_ribbon(alpha =0.6, color =NA, fill ="orange") +geom_line() +labs(x ="FWI", y =expression(lambda)) +# Agregamos los datosgeom_point(aes(fwi, fires), data = datos, inherit.aes = F,alpha =0.7, size =2)
Code
# La cinta más ancha es donde esperamos ver los datos con 95 % de probabilidad;# la cinta más fina es donde esperamos ver la media de muchos datos, con 95 % # de probabilidad.# El más ancho se suele llamar intervalo de predicción; el más fino, intervalo de # credibilidad para la media.
\(R^2\) bayesiano (Gelman et al. 2019)
El coeficiente de determinación bayesiano se define así:
\[
R^2 = \frac{\text{varianza de los predichos}}{\text{varianza de los predichos} + \text{varianza residual}},
\]
donde los predichos son la esperanza (media) predicha para cada observación. La varianza de los predichos se toma a lo largo de las observaciones (\(N\) predichos). La varianza residual es la varianza estimada para cada observación. La varianza residual se promedia entre las \(N\) observaciones.
Representa la proporción de la variabilidad en la distribución estimada de los datos que es explicada por el modelo.
\(R^2\) a partir de muestras de la posterior
Con una estimación bayesiana basada en muestras de la posterior \(s \in \{1, \dots, S\}\), el \(R^2\) se puede calcular para cada muestra:
donde \(\boldsymbol{\mu}_s\) es el vector de largo \(N\) con la esperanza estimada para cada observación, calculada con el vector de parámetros \(\boldsymbol{\theta}_s\). \(\boldsymbol{\tau}_s\) es el vector de varianzas estimadas para cada observación, también basado en \(\boldsymbol{\theta}_s\). En el caso de un modelo normal, sería un vector de \(\sigma^2\). (Nótese que en la mayor parte de las distribuciones la varianza varía con la media, por lo que suele variar entre observaciones.)
Estimando \(R_s^2\) para cada \(s\), obtenemos su distribución posterior.
Importante: dado que el \(R^2\) bayesiano está basado en la varianza y media estimadas, para que sea confiable, el modelo debe representar los datos de forma adecuada. Por ejemplo, un modelo que subestima la variabilidad en los datos va a sobreestimar el \(R^2\). Por eso es importante verificar el modelo con residuos DHARMa y hacer lo posible por que cuantifique de forma razonable la variabilidad.
\(R^2\) para el modelo de incendios.
Code
# Vector para almacenar el R2R2 <-numeric(S)# Loop sobre muestras de la posteriorfor (i in1:S) {# predichos para cada observación mu <-exp(d$a[i] + d$b[i] * datos$fwi) # es lambda!# varianza de los predichos = varianza explicada var_fit <-var(mu)# varianza residual para cada observación. Esto depende de la # distribución. En la poisson, la varianza es igual al media var_res <- mu # sólo para ser explícitos# y tomamos su media entre observaciones var_res_mean <-mean(var_res)# R2 R2[i] <- var_fit / (var_fit + var_res_mean)}plot(density(R2, from =0, to =1), main =NA, xlab ="R2")
El \(R^2\) da alto, pero hay que considerar que el modelo subestima la variabilidad en los datos (variabilidad no explicada). Si modeláramos con una binomial negativa, que admite más variabilidad que una poisson, el \(R^2\) sería notablemente más bajo.
Formulación del \(R^2\) bayesiano:
Gelman, A., Goodrich, B., Gabry, J., & Vehtari, A. (2019). R-squared for Bayesian regression models. The American Statistician.
Processing csv files: /tmp/RtmpVB3dXY/germinadas-202505052113-1-607bac.csv, /tmp/RtmpVB3dXY/germinadas-202505052113-2-607bac.csv, /tmp/RtmpVB3dXY/germinadas-202505052113-3-607bac.csv, /tmp/RtmpVB3dXY/germinadas-202505052113-4-607bac.csv
Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.
Checking sampler transitions for divergences.
No divergent transitions found.
Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.
Effective sample size satisfactory.
Split R-hat values satisfactory all parameters.
Processing complete, no problems detected.
Verificación predictiva posterior
Code
# Extraemos muestras de los parámetros que necesitamosd <- fit$draws(c("alpha", "beta_D", "beta_E", "beta_H", "phi"), format ="draws_df")colnames(d) <-c("a", "bD", "bE", "bH", "phi")S <-nrow(d)N <-nrow(datos)# Verificamos el ajuste simulando de la predictiva posterior. Para simular de # la beta-binomial en R usaremos el paquete extraDistr, que incluye la función# rbbinomysim <-matrix(NA, N, S)for (i in1:S) { mu <-plogis( # equivale a inv_logit = logit inverso d$a[i] + d$bD[i] * datos$D + d$bE[i] * datos$E + d$bH[i] * datos$H ) a <- mu * d$phi[i] b <- (1- mu) * d$phi[i] ysim[, i] <-rbbinom(N, datos$sembradas, a, b)} ysim_t <-t(ysim)# Miramos con bayesplotppc_intervals( datos$germinadas, yrep = ysim_t, x = datos$danio,prob =0.5, prob_outer =0.9)
Code
use <-sample(1:S, 100)ppc_dens_overlay(datos$germinadas, yrep = ysim_t[use, ])
plotResiduals(res, form = datos$altitud, rank = F)
Code
plotResiduals(res, form = datos$altura, rank = F)
Comparamos posterior vs. previa de \(\phi\)
Code
# phiplot(density(d$phi, from =0), main =NA, xlab =expression(phi))curve(dlnorm(x, log(7), 1.5), add = T, col =4) # previa definida en Stan
\(R^2\)
Code
# Vector para almacenar el R2R2 <-numeric(S)# Por comodidad, definimos las sembradasSS <- datos$sembradas[1]# Loop sobre muestras de la posteriorfor (i in1:S) {# En la beta-binomial, la media es mu * size,# siendo size la cantidad de sembradas# predichos para cada observación mu <-plogis( # equivale a inv_logit = logit inverso d$a[i] + d$bD[i] * datos$D + d$bE[i] * datos$E + d$bH[i] * datos$H ) media <- mu * SS # 120 sembradas# varianza de los predichos = varianza explicada var_fit <-var(media)# varianza residual para cada observación. En la beta-binomial, la # varianza se define así: var_res <- SS * mu * (1- mu) * ((N + d$phi[i]) / (1+ d$phi[i]))# y tomamos su media entre observaciones var_res_mean <-mean(var_res)# R2 R2[i] <- var_fit / (var_fit + var_res_mean)}plot(density(R2, from =0, to =1), main =NA, xlab ="R2")
Visualizamos resultados de interés: predicciones parciales
Code
# Predicción parcial para visualizar el efecto del daño por fuego (D).# Graficaremos curvas basadas en nrep valores de daño en el rango observadonrep <-200danio_seq <-seq(min(datos$danio), max(datos$danio), length.out = nrep)# Estandarizamos la secuencia, porque los parámetros viven en escala estandarizada.D_seq <- (danio_seq - D_mean) / D_sd# Creamos data.frame de predicciones con las predictoras. Las predictoras no # focales se fijan en su media (0 porque están estandarizadas), por eso es una # predicción parcial.pred <-data.frame(D = D_seq, E =0,H =0)pred$danio <- danio_seq # la agregamos en la escala original, para visualizar# Matriz para llenar con predicciones de mumu_mat <-matrix(NA, nrep, S)# Y matriz para llenar con valores simulados de semillas germinadasy_mat <-matrix(NA, nrep, S) # se parece a ysim, pero tiene nrep filas, no N.# Las llenamosfor (i in1:S) { mu <-plogis( # equivale a inv_logit = logit inverso d$a[i] + d$bD[i] * pred$D + d$bE[i] * pred$E + d$bH[i] * pred$H # los dos últimos términos pueden obviarse, pero los dejamos para # explicitar. ) mu_mat[, i] <- mu# calculamos alpha y beta para simular a <- mu * d$phi[i] b <- (1- mu) * d$phi[i] y_mat[, i] <-rbbinom(nrep, datos$sembradas[1], a, b)} # Resumimos la distribución posterior de mu fila por filapred$mu_mean <-rowMeans(mu_mat)pred$mu_lwr <-apply(mu_mat, 1, quantile, probs =0.025)pred$mu_upr <-apply(mu_mat, 1, quantile, probs =0.975)# Convertimos la distribución de y en proporciones, para que esté en la misma# escala que muy_mat_p <- y_mat / datos$sembradas[1]# y resumimospred$y_lwr <-apply(y_mat_p, 1, quantile, probs =0.025)pred$y_upr <-apply(y_mat_p, 1, quantile, probs =0.975)# Calculamos la proporción observada en los datosdatos$prop <- datos$germinadas / datos$sembradas# Graficamos ggplot(pred, aes(danio, mu_mean, ymin = mu_lwr, ymax = mu_upr)) +# Intervalo de predicción (posterior predictive)geom_ribbon(aes(danio, ymin = y_lwr, ymax = y_upr), inherit.aes = F, color =NA, alpha =0.3, fill ="blue") +# Intervalo para la mediageom_ribbon(color =NA, alpha =0.7, fill ="orange") +# Media de la mediageom_line() +# datosgeom_point(aes(danio, prop), data = datos, inherit.aes = F,size =2, alpha =0.7) +labs(y ="Proporción de germinación",x ="Daño (%)" )
En esta predicción se fijaron en su media las demás predictoras. Esto hace que la predicción no sea comparable con los datos, ya que en los datos, las predictoras sí varían. En este caso, los datos no discrepan de la predicción porque las otras predictoras tienen efectos muy pequeños.
Para realizar otras predicciones hay que ponerse creativo al crear la tabla con nuevos valores para las predictoras. Por ejemplo, podríamos graficar más de una curva, para 3 valores contrastantes de alguna otra predictora. Esto es útil para visualizar interacciones.
Simulaciones previas
Al definir las previas de un modelo con ajuste bayesiano necesitamos inicialmente posicionarnos en relación a qué efecto deseamos de las previas. Algunas opciones pueden ser:
reflejar conocimiento previo sobre el sistema,
restringir levemente los parámetros para que tomen valores razonables,
restringir fuertemente los parámetros para que el modelo sea estimable,
influir lo menos posible en el resultado, limitándose a garantizar la convergencia de algoritmos de estimación.
Independientemente de la postura, para definir previas necesitamos poder traducir información de formas variadas a una distribución de probabilidad. Salvo que contemos con una intuición algebraica y probabilística envidiable, necesitamos graficar las implicancias de las previas en los modelos. Claramente, en contextos donde no hay información previa abundante, conviene investigar un poco qué tipo de previas se recomiendan para el tipo de modelo en cuestión. Para muchas estructuras hay recomendaciones disponibles.
En última instancia, si la definición de previas nos quita el sueño, podemos ajustar el modelo con distintas previas y comparar resultados. Esto incluso puede plantearse como un problema de comparación de modelos: cambiar la previa es equivalente a cambiar el modelo.
Para mostrar estrategias para definir previas utilizaremos los modelos ajustados arriba para el número de incendios y la cantidad de semillas germinadas. Supondremos que no hay mucha información previa y que deseamos distribuciones previas con poco efecto en el ajuste, de tal forma que la posterior esté dominada por la verosimilitud.
Si deseamos previas que generen valores razonables de \(\lambda\) sin afectar mucho el resultado, estamos buscando previas levemente regularizadoras.
Un buen primer paso es intentar interpretar los parámetros del modelo. Si eso es muy difícil, conviene graficar.
Consideremos que \(\alpha = \log(\lambda)\) cuando \(FWI = 0\). De esta manera, \[\lambda_{FWI=0} \sim \text{log-normal}(\mu_{\alpha}, \sigma_{\alpha}).\] Entonces, graficamos la lognormal variando sus parámetros (expresados en escala log).
Para graficar previas o sus implicancias rápidamente, resulta de gran utilidad la función curve, que evalúa una expresión conteniendo “x”.
Code
mu_a <-log(10) # la mediana de la lognormal sería 10sigma_a <-2curve(dlnorm(x, mu_a, sigma_a), from =0, to =50)
Con \(\sigma_{\alpha} = 2\) la distribución se hace muy amplia. No nos dejemos engañar por su pico cerca de cero. Si bien en esa zona la densidad es alta, la masa de probabilidad (área bajo la curva) es baja, por lo que esta previa no sesgaría la estimación hacia cero si hay al menos unos pocos datos.
Pero si nos interesa no sesgar el resultado, otra estrategia es centrar la previa en la media de los datos, y darle un desvío grande:
Code
# Cargamos datosdatos <-read.csv(here::here("R_curso_modelos", "datos","barbera_data_fire_total_climate.csv"))mu_a <-log(mean(datos$fires)) sigma_a <-2curve(dlnorm(x, mu_a, sigma_a), from =0, to =50)
Si queremos informar aún menos, podemos aumentar un poco más $_{}.
Ahora nos centramos en \(\beta\). Por suerte, para este modelo (GLM con enlace log), \(\beta\) tiene una interpretación clara: \(\exp(\beta)\) es el factor en que aumenta \(\lambda\) cuando \(x = FWI\) aumenta una unidad. Si suponemos que \(\lambda\) se duplica al aumentar el \(FWI\) una unidad, \(\beta = \log(2)\). Pero como el \(FWI\) varía entre 0 y ~50, ese aumento sería absurdamente grande. Quizás es más sensato pensar en aumentos más pequeños. Probemos un aumento de 1.1 \(\beta = \log(1.1)\).
Nuevamente, \(\exp(\beta)\) sigue una distribución log-normal:
Code
mu_b <-log(1.1)sigma_b <-1curve(dlnorm(x, mu_b, sigma_b), from =0, to =5)
Esta previa tiene mediana y media positiva (1.1), pero también permite valores menores a 1, que implican disminución del número de incendios con el \(FWI\). Si esto no nos parece razonable, podemos modificar la previa para dejar poca probabilidad abajo de 1 (restricción blanda) o truncar la previa en 1, para que sólo permita aumentos. Iremos por el primer enfoque.
Code
mu_b <-log(1.1)sigma_b <-0.05curve(dlnorm(x, mu_b, sigma_b), from =0, to =5)
Ahora es una previa más razonable, quizás algo informativa. Pero cuán informativa es dependerá de cuán informativos sean los datos.
Ahora graficamos las previas elegidas para \(\alpha\) y \(\beta\) (escala log, son normales):
Code
par(mfrow =c(1, 2))curve(dnorm(x, mu_a, sigma_a), from =-10, to =15, xlab =expression(alpha),ylab ="Densidad")curve(dnorm(x, mu_b, sigma_b), from =-1, to =1, xlab =expression(beta),ylab =NA)
Code
par(mfrow =c(1, 1))
Ahora podemos simular valores de estas previas para explorar qué curvas de \(\lambda = f(FWI)\) implican:
Code
light_col <-rgb(0, 0, 0, 0.1)# gráfico iniciala <-rnorm(1, mu_a, sigma_a)b <-rnorm(1, mu_b, sigma_b)curve(exp(a + b * x), from =0, to =20, ylim =c(0, 100), col = light_col,xlab ="FWI", ylab =expression(lambda)) # Repetimos con muchas muestras, usando el argumento "add = TRUE" en curve()for (i in1:200) { a <-rnorm(1, mu_a, sigma_a) b <-rnorm(1, mu_b, sigma_b)curve(exp(a + b * x), col = light_col, add = T)}
Muchas curvas parecen poco razonables. Pero bueno, estábamos buscando previas poco informativas, ¿no?
También podemos simular datos utilizando los valores observados de la predictora. Esto implica simular de la distribución predictiva previa (análoga a la distribución predictiva posterior).
Code
N <-nrow(datos)S <-10000ysim <-matrix(NA, N, S)for (i in1:S) { a <-rnorm(1, mu_a, sigma_a) b <-rnorm(1, mu_b, sigma_b) lambda <-exp(a + b * datos$fwi) ysim[, i] <-rpois(N, lambda)}hist(as.numeric(ysim))
Los datos simulados pueden tomar valores muy altos. Esto puede ser razonable si no queremos que la previa influya mucho en los resultados.
Para ilustrar un ejemplo de previa poco intuitiva, alteraremos el modelo: en vez de suponer una distribución poisson, supondremos binomial negativa, que admite mayor dispersión mediante el parámetro \(phi\):
Por lo que \(\phi\) tiene una relación inversa con la varianza, y se puede decir que es un parámetro de precisión (opuesto a dispersión).
Podemos evaluar las implicancias de \(\phi\) fijando \(\lambda\) y visualizando la pmf.
Code
# usamos dnbinom de extraDistr, que la parametriza igual que Stan's neg_binomial_2.# size = phix <-0:100y <-dnbinom(x, mu =10, size =1000)barplot(y ~ x); abline(v =10, col ="red", lty =2, lwd =1.5)
Code
# Con size alto tiende a poissony <-dnbinom(x, mu =10, size =1)barplot(y ~ x); abline(v =10, col ="red", lty =2, lwd =1.5)
Code
y <-dnbinom(x, mu =10, size =0.1)barplot(y ~ x); abline(v =10, col ="red", lty =2, lwd =1.5)
Code
# ChatGPT recomendó una previa gama(2, 0.1), lo que implica# a = 2 = 1 / 0.5# b = 0.1 = 1 / (0.5 * 20) # lo que implica mu_phi = 20, tau = 0.5lambda <-20# probar cambiarlomu_phi <-20# nos gustan?tau <-0.5phi <-rgamma(1, 1/ tau, 1/ (mu_phi * tau))curve(dnbinom(x, size = phi, mu =10), # este mu es lambdato = lambda *10, col = light_col, ylim =c(0, 0.5)) for (i in1:200) { phi <-rgamma(1, 1/ tau, 1/ (mu_phi * tau))curve(dnbinom(x, size = phi, mu =10), add = T, col = light_col) # este mu es lambda}
Code
# Para ser menos informativo, prefiero mu_phi = 20, tau = 1 (tau es la dispersión# de la gama, por lo que mi previa es más amplia que la recomendada por ChatGPT)
Simulaciones previas 2: germinación
Retomemos el modelo de semillas germinadas de Alinari et al. 2024:
Con predictoras estandarizadas (media = 0), \(\alpha\) es la probabilidad de germinación (\(\mu\)) cuando las predictoras están en sus medias.
Si \(\alpha \sim \text{normal}()\), entonces \(\text{logit}^{-1}(\alpha) \sim \text{logit-normal}()\).
Si la media de la previa de \(\alpha\) es 0, a medida que aumentamos su desvío acumulamos masa cerca de \(\mu = 0\) y \(\mu = 1\):
Code
# usamos el paquete logitnorm para la distribución logit-normalmu_a <-0sigma_a <-1.4curve(dlogitnorm(x, mu_a, sigma_a))
Code
sigma_a <-1.8curve(dlogitnorm(x, mu_a, sigma_a))
Code
sigma_a <-10curve(dlogitnorm(x, mu_a, sigma_a))
Si bien puede resultar poco intuitivo, la previa con sd = 10 resulta muy poco informativa. Permite que se estimen probabiliades extremas si los datos lo indican.
Yendo a los \(\beta\), una gran ventaja de estandarizar las predictoras es que sus previas pueden ser todas iguales, a falta de información que sugiera que lo sean.
Para explorar previas sobre \(\beta\) conviene asumir el efecto parcial de una predictora, fijando las demás en sus medias y \(\alpha = 0\).
Como seremos indiferentes al signo, las previas de \(\beta\) estarán centradas en 0.
Code
sigma_b <-0.5b <-rnorm(1, 0, sigma_b)# plogis es inv_logit en R.# de -3 a 3 porque ese rango suele ocupar una variable estandarizada# cuando es aproximadamente normal.curve(plogis(b * x), from =-3, to =3, ylim =c(0, 1),col = light_col)for (i in1:200) { b <-rnorm(1, 0, sigma_b)curve(plogis(b * x), col = light_col, add = T)}
Code
# sigma mayorsigma_b <-2b <-rnorm(1, 0, sigma_b)curve(plogis(b * x), from =-3, to =3, ylim =c(0, 1),col = light_col)for (i in1:200) { b <-rnorm(1, 0, sigma_b)curve(plogis(b * x), col = light_col, add = T)}
Code
# y mayorsigma_b <-5b <-rnorm(1, 0, sigma_b)curve(plogis(b * x), from =-3, to =3, ylim =c(0, 1),col = light_col)for (i in1:200) { b <-rnorm(1, 0, sigma_b)curve(plogis(b * x), col = light_col, add = T)}
Igual que con \(\alpha\), una previa amplia para \(\beta\) no sesga la posterior hacia valores extremos; simplemente no impide que los datos lo sugieran. No es conveniente si tenemos escasos datos.
Por último, debemos definir la previa para \(\phi\), el parámetro de precisión de la distribución beta-binomial. La beta-binomial es una binomial con mayor dispersión (regulada por \(\phi\)), que surge de suponer que la probabilidad de éxito (\(\pi\), aquí \(\mu\)) de una binomial sigue una distribución beta. Su varianza es
donde \(N\) es el número de experimentos, \(\phi\) es un parámetro de precisión (\(\alpha + \beta\) en la distribución beta) y \(\mu\) es la probabilidad de éxito. (La media es \(\mu \times N\).)
Graficaremos la pmf suponiendo \(N = 500\) y \(\mu = 0.5\).
Code
# usamos extraDistr# dbetabinom, se parametriza con alpha y beta, no con mu y phi.# Se definen así:mu <-0.5N <-500# Nphi <-rlnorm(1, log(10), 2)alpha <- mu * phi # necesarios para dbbinombeta <- (1- mu) * phicurve(dbbinom(x, size = N, alpha = alpha, beta = beta),to = N, col = light_col, ylim =c(0, 0.5)) for (i in1:200) { phi <-rlnorm(1, log(10), 2) alpha <- mu * phi beta <- (1- mu) * phicurve(dbbinom(x, size = N, alpha = alpha, beta = beta),col = light_col, add = T)}
La previa definida es muy amplia; permite muchas distribuciones en que los valores se acumulan en los extremos (todos éxitos o todos fracasos). Podemos encogerla.
Code
# usamos sd = 1.5 en vez de 2 en la previa lognormal de phiphi <-rlnorm(1, log(10), 1.5)alpha <- mu * phi # necesarios para dbbinombeta <- (1- mu) * phicurve(dbbinom(x, size = N, alpha = alpha, beta = beta),to = N, col = light_col, ylim =c(0, 0.1)) for (i in1:200) { phi <-rlnorm(1, log(10), 1) alpha <- mu * phi beta <- (1- mu) * phicurve(dbbinom(x, size = N, alpha = alpha, beta = beta),col = light_col, add = T)}
